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Abstract 

Motivated by indications that heavy (charm and bottom) quarks interact strongly at temper- 
atures generated in heavy ion collision experiments, we suggest a non-perturbative definition 
of a heavy quark chemical equilibration rate as a transport coefficient. Within leading-order 
perturbation theory (corresponding to 3- loop level), the definition is argued to reduce to 
an expression obtained from the Boltzmann equation. Around T ~ 400 MeV, an order-of- 
magnitude estimate for charm yields a rate rT > 60 fm/c which remains too slow to play 
a practical role in current experiments. However, the rate increases rapidly with T and, due 
to non-linear effects, also if the initial state contains an overabundance of heavy quarks. 
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1. Introduction 



In a fully thermalized medium, the momenta of bosons and fermions are distributed accord- 
ing to the Bose and Fermi distributions, respectively, parametrized by a single temperature, 
T, and chemical potentials associated with conserved global charges. In contrast, the most 
important cosmological relics, such as Light Element Abundances, Dark Matter, or Baryon 
Asymmetry, rely on deviations from thermal equilibrium. In a canonical Dark Matter sce- 
nario, for instance, the overall abundance of the Dark Matter particles is determined through 
a "freeze-out" period, which takes place when their annihilation rate becomes too slow to 
track the total number density determined by the Fermi distribution, which decreases expo- 
nentially when irT <C M, where M denotes the particle mass. Since the number densities of 
particles and antiparticles remain equal, this deviation cannot in relativistic field theory be 
represented through a chemical potential, and we speak of chemical non-equilibrium. (Typ- 
ically, elastic scatterings with the plasma particles still continue after this period, so that 
kinetic equilibrium is maintained down to lower temperatures, cf. e.g. ref. [1].) A freeze-out 
process leading to chemical non-equilibrium is also responsible for the ~ 20% primordial 
helium abundance observed in the Universe today, cf. e.g. ref. [2]. 

Analogous processes are assumed to play a role in heavy ion collisions. In particular, for 
7rT <C M, the kinetic equilibration rate of heavy quarks scales as T^m ~ a 2 s \n{a s )T 2 / M [3]- 
[6], whereas the chemical equilibration rate scales as r c h em ~ a^Ts exp(— M/T)/Mz [8, 9]. 
Experimental data from RHIC and LHC suggest that charm quarks do have time to kinetically 
equilibrate, thereby participating in hydrodynamic flow (cf. e.g. refs. [10, 11]), and theoretical 
efforts to understand this up to the non-perturbative level are under way [12]-[14]. Building 
on earlier studies of strange quarks [15] it is believed, in contrast, that chemical equilibration 
does not take place; the number density of charm quarks and antiquarks is essentially assumed 
to remain as determined by an initial hard process [16], implying that there are more heavy 
quarks present than would be due for chemical equilibrium (cf. e.g. ref. [17]). 

The purpose of this study is to suggest a definition of a chemical equilibration rate of heavy 
quarks near equilibrium, similarly to what was achieved for their kinetic equilibration rate 
earlier on [18, 19]. A definition should be possible in the heavy-quark limit M S> vrT, in 
which the rate itself is much slower than typical "fast" plasma rates, Ff ast ~ a™T, n > 1. 
(If no scale separation is present between M and ttT, then pair creations and annihilations 
take place as fast as elastic processes, and the massive degrees of freedom are to a good 
approximation in full thermal equilibrium with the strongly interacting heat bath.) 

The plan of this paper is the following. After some general considerations in sec. 2, we 
recall the derivation of the chemical equilibration rate to leading order in a s , making use of 
the Boltzmann equation, in sec. 3. This is followed by a reminder that loop corrections are 
likely to be substantial at any realistic temperature, in sec. 4. A non-perturbative formulation 
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is put forward in sec. 5. Subsequently we argue, in sec. 6, that in the weak-coupling limit the 
expression of sec. 5 reduces to the result of sec. 3. A brief discussion of implications as well 
as prospects for non-perturbative studies concludes this writeup in sec. 7. 

2. General considerations 

Assume that the system possess an approximately conserved particle number. Let us denote 
the corresponding number density 1 by n(t). In thermal equilibrium the value of n fluctuates 
around its equilibrium value. To treat the non-equilibrium problem we follow the general 
method described in ref. [7]. Let 5n(t) = n(t) — n eq at some time t be large compared to the 
mean fluctuation. It will then evolve towards its equilibrium value. Let us assume that the 
characteristic time scale r for this evolution is much larger than the other relaxation times of 
the system. We only want to resolve time scales of order r. Then the non-equilibrium state 
is completely characterized by the instantaneous value of Sn. Therefore the time derivative 
of Sn can only depend on the value of Sn and on thermodynamic quantities of the system 
such as temperature and chemical potentials. When 5n is sufficiently small, one can expand 
Sn in powers of Sn and keep only the linear term, 

Sn(t) = -T chcm Sn(t) . (2.1) 

The coefficient r chcm only depends on thermodynamic quantities. 

Let us now be specific and choose n to be the sum of quark and antiquark number densities, 

n = n Q + n-Q . (2.2) 

We consider the heavy quark baryon number density n Q — n^ to vanish (i.e. the baryon 
chemical potential to be zero). We are interested in the limit that irT <C M. For heavy 
particles, {Sn(t)}\ oss ~ e -2M / T , because a heavy quark-antiquark pair gets annihilated, and 
Sn(t) ~ n cq ~ e~ M l T . Therefore T itself scales as ~ e~ M / T , implying that this rate is much 
slower than most other processes in the system. In particular, this rate is slower than the 
kinetic equilibration rate. Therefore the heavy quarks can be considered to be in kinetic 
equilibrium, which means that they move very slowly. These almost static quarks expe- 
rience rare number changing reactions, and a non-perturbative description of the resulting 
dynamics, incorporating both the non-equilibrium evolution of eq. (2.1) as well as equilibrium 
fluctuations, is presented in eqs. (5.11)-(5.20) below. 

x It is important to consider the number density rather than the differential phase space distribution, 
because otherwise it would be difficult to distinguish between processes changing the kinetic and the chemical 
decomposition of the system. 
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3. Boltzmann equation 



If the system is weakly coupled, one can usually compute the coefficient r chem in eq. (2.1), at 
least to leading order, from the Boltzmann equation. If we take into account 2—7-2 scattering 
processes and consider the limit 7rT <C M, it takes the form (cf. e.g. ref. [20]) 



h = -c(n - n cq ) = n loss + n gain , (3.1) 

where ni oss = —cn 2 . In equilibrium, with n(t) = n cq , gain and loss terms must cancel each 
other, and the number density is constant. Now linearize (3.1) as described in sec. 2, which 
gives 5ii = —2cnSn. Thus we can obtain r chcm from the loss term in eq. (3.1) via 

r chem = -2^. (3.2) 

An analogous discussion, implemented by introducing separate "chemical potentials" for the 
quarks and antiquarks, can be found in ref. [15]. 

Now we compute r chem using eq. (3.2) with tree-level matrix elements. The relevant loss 
processes are shown in fig. 1. Inserting the number of degrees of freedom of the initial state, 
2N C , the decay rate according to eq. (3.2) can be written as 

x fAEk 1 )fAE ka )l\T l \M 1 \ 2 [l + f B {e pi )] [1 + /bM] 

+ N f j:\M 2 \ 2 [l-f F (e Pl )] [1 - / F (e P2 )] } . (3.3) 

Here J k = J ^p-; k a are momenta in the initial state and pi those in the final state; E^ a = 
\/k 2 + M 2 is the energy of a massive particle and e Pi = |pj| is that of a massless one; 
and / f ,/b are the Fermi and Bose distributions, respectively. The sums are taken over the 
quantum numbers of all on-shell degrees of freedom, i.e. 2N C for quarks and antiquarks, and 
2d a for gluons, with cLa = N 2 — 1. By iVf we denote the number of light quark flavours, and 
later on Cf = cIa/(2N c ) will also appear. The factor | in front of the gluonic amplitude 
accounts for the two final state particles being identical [15]. 

Taking the amplitude M. 2 °f fig- 1 as an example, a text-book calculation yields (cf. e.g. 
refs. [21, 22]) 

£ |7W 2 | 2 = ^ 4 CfN c r (M2 _ 2 + 2 _ 2 + 2M 2 1 j (3 4) 

where s,t,u are the standard kinematic invariants: s = (V\ + V2) 2 = (/Ci + /C2) 2 ; t = 
(Pi - /d) 2 = (P 2 - /C 2 ) 2 ; and u = (Pi - /C 2 ) 2 = (P 2 - /d) 2 . 
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Figure 1: Scatterings through which an overabundance of heavy quarks can disappear, assuming that 
there is an exponentially small thermal distribution of antiquarks present (or vica versa). A double 
line indicates heavy quarks, a single line light quarks, and a wiggly line gluons. 



The result simplifies further in the heavy-quark limit. Because of Boltzmann suppression 
of f F (Ek a ) at M S> 7rT, we can consider the decaying heavy quark and antiquark to be almost 
at rest with respect to the thermal medium: 



/C 2 



k 2 

M + -^_,k 2 
2M' 



(3.5) 



with k a ~ V ttTM <C M. In contrast p\ and P2 are large because they have to carry away 
the energy liberated in the pair annihilation. So ki + k 2 can be approximated as zero in the 
phase space constraints, and the Fermi distributions fw( e pi) can be omitted: 

e -M/T , d 3 k 



chem 



2 fc 2 
e 2MT 



4iV c M 2 J (2vr) 3 
1_ /-d 3 pi /-d 3 P2 

(27TV 



5(3)(pi + P2)5(epi + £pa " 2M) ^ |Als 



(3.6) 



Here we cancelled a factorized integral against the one in the denominator. Noting also that 
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-M 



-M' L 



(3.7) 



we get ^2 \-M-2 1 2 ~ ^9 A CfN c . The remaining integrals are trivially carried out, and we obtain 



{qq) ^ g CfN{ / TM \ f m/t 

chem ~ 8yrM 2 \ 2 TT ) 



(3.8) 



A similar computation can be carried out with gluons, represented by the amplitude Ai 1 
of fig. 1. Again the result is well-known (cf. e.g. refs. [21, 22]), and reads 



El-Mil 



4g*C F N c < 4N C 
+ 2C F 
- 2N r 



(M 2 -t)(M 2 



(M 2 -t)(M 2 



+ (2C F - N c 



2M 2 (s-AM 2 ) 
(M 2 - t)(M 2 - u) 



2M 2 (M 2 + t) 



{M 2 - t) 2 
(M 2 - t)(M 2 - it) + M 2 {u - t) 



+ (t4->u) 



s(M 2 - t) 



(3.9) 
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Figure 2: Examples of 1-loop corrections to the scattering amplitude M. 2 °f n §- 1- 



In the heavy-quark limit, eq. (3.7), this simplifies to £ {M^ 2 ^ 4g 4 C F N c (4:C F - N c ). The 
phase space integration goes through as before, and recalling the | in eq. (3.3), eq. (3.8) gets 
completed into 



chcm 



8ttM 2 V 1 ' 2 A 2tt 7 



-M/T 



(3.10) 



Numerically 2C F - iVc/2 = \ for 7V C = 3; for N { = this agrees with eq. (10) of ref. [9]. (We 
note, however, that for three light flavours, i.e. Nf = 3, fermionic final states are significantly 
more important than purely gluonic ones.) 



4. Towards loop corrections 

The result of eq. (3.10) could well suffer from large radiative corrections. A few representative 
examples of next-to-leading order (NLO) amplitudes are shown in fig. 2. In particular, the 
first amplitude, iterated by further rungs connecting the heavy quark and antiquark to each 
other, is responsible for binding the particles to a quarkonium-like resonance. In the context 
of Dark Matter co-annihilation, such a threshold enhancement is assumed to play a potentially 
important role, cf. e.g. refs. [23, 24]. However, this is not the only class of processes in our 
case: as illustrated in fig. 2, all participating particles carry a colour charge, so that there 
may also be final-state interactions, as well as "non-factorizable" terms connecting the initial 
and final states. 

For future reference, we remark that there is one Euclidean observable in which rungs 
between the heavy particles can also appear but which is nevertheless very well understood. 
This is the heavy quark-number susceptibility, formally defined as 

X/ = y((^7o^)(r,x)(^7oV)(0,0)) T , < r < (3 , = i . (4.1) 

Because of charge conservation the argument r can be chosen at will. With vanishing chemical 
potentials, the susceptibility measures the mean number of heavy particles created by thermal 
fluctuations, and is therefore closely related to the distribution function f F (Ek 2 ) on which 
the heavy quarks scatter in eq. (3.3). 

We recall that in the free limit the susceptibility evaluates to 

Xf = m -Jj^y f*(Ek) [1 " fAEk)] ■ (4-2) 
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For massless quarks the integral can be carried out in a closed form, yielding Xf = N c T 3 /3, 
to which loop corrections are known up to a high order [25], generically decreasing the sus- 
ceptibility from the free value. To us more relevant is the non-relativistic limit, 




Here the temperature dependence is precisely the same as that in eq. (3.8). Lattice data 
indicate that the susceptibility grows rapidly with the temperature and, in the charm case, 
overcomes the exponential suppression already at temperatures of a few hundred MeV [26]- 
[28] , in line with the general expectation [29] . We will keep these observations in mind when 
estimating the numerical importance of the exponential suppression in sec. 7. 

5. Non-perturbative formulation 

Motivated by the remarks in sec. 4, the goal now is to suggest a non-perturbative definition of 
the heavy quark chemical equilibration rate. This could allow for a systematic computation 
of higher order corrections, or in principle be subjected e.g. to a lattice investigation. 

In relativistic theories there is no obvious definition for a particle number operator. Here 
we are interested in heavy quarks and antiquarks with very small velocities. In this case 
the energy of quarks and antiquarks is roughly given by the sum of their rest energies or, in 
other words, by their number density times the heavy quark mass M. Therefore the energy 
density of heavy quarks and antiquarks is a good measure for their number density. We 
propose to define the relaxation time of the number density n = n Q + n-^ through the real 
time correlation function of the heavy quark Hamilton operator. 

We start by introducing an operator describing heavy quark energy loss, both through 
elastic and through inelastic processes (sec. 5.1); define then a "transport coefficient" related 
to this operator, capturing the desired rate (sec. 5.2); and finally simplify one of the correlators 
appearing by considering the heavy-quark limit (sec. 5.3). 

5.1. Operator for heavy quark energy loss 

A form of the fermionic energy-momentum tensor which is symmetric, gauge-invariant, and 
leads to a correct finite trace anomaly, reads [30, 31] 

2f = ^^f^ + Y^^-rT^- (5-1) 
Here rf v = diag(H ) and 

^y^V = ^y^V - ^ylJ^v > ( 5 - 2 ) 
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with L^ u ip = (d u — igA v )tp, tp 1)^ = tp(d u + igA u ), and g denoting the bare gauge coupling. 
The Lagrangian can be written with a similar notation as 



(5.3) 



The heavy quark Hamilton operator is now defined by taking a spatial integral over T f 00 , 
with the fields promoted to operators: 



H 



55 /x ff00 = jj{- l 2^ + M )^- 



(5.4) 



Summation over repeated spatial indices is understood. Obviously, H could be written in 
other forms by use of the Dirac equation, but for us it appears to be beneficial to employ a 
version with spatial derivatives only, because then partial integrations are formally allowed. 

In order to derive the operator for energy loss, let us also write down the Dirac equation 
in an explicit form, by placing time derivatives on the left-hand side: 



d t i> 



-i(M 7 ° - gA ) - j°7 j Dj 



4> 



i(M 7 ° - gA ) 



(5.5) 
(5.6) 



In all of what follows, equations of motion are used for fermions only; derivatives acting on 
gauge fields are left "as is" , formally assuming that gauge fields form a differentiable off-shell 
background over which a path integral is to be carried out at a later stage. 

The task now is to construct dfH. The derivative can act on any of the three possible 
locations in eq. (5.4): 



d t H = J I (dt$) {-i^Dj +M^ + 



4, -gj j d Aj U + M i-/ j D + M ) ( d t ip 



(5.7) 



Inserting eqs. (5.5), (5.6) and carrying out one partial integration, numerous cancellations 
take place, and we are finally left with 



8 t H = -9 J ^(doAj-djAo-igAoAj + igAjAo^ = -g J 



g / V I 3 F 0j i> 



(5.8) 



So, in the presence of interactions (g 7^ 0), the energy carried by heavy quarks is not conserved. 

It appears that eq. (5.8) has a classical interpretation. If a charged particle feels a Lorentz 
force, 

^ = ?(e + vxb), (5.9) 

then its energy changes as 



dE 



dp 



V ^'d* 



d P 

v- = ,v.E. 



(5.10) 
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Recalling that tpj^ are the spatial components of a current this is seen to agree in form 
with eq. (5.8). However, being a Fock space operator, dfH of eq. (5.8) describes also number- 
changing reactions; in particular, if the initial state has more quarks and antiquarks than 
would be due for chemical equilibrium, a net pair annihilation should take place, and in the 
large-time limit the corresponding matrix elements dominate the statistical average of dtH. 

5.2. Defining a transport coefficient 

To describe the depletion of an overabundance of heavy quarks through a single coefficient, 
we follow a general method which has also been used for determining their kinetic equili- 
bration rate [18, 19]. The goal is to relate the non-equilibrium rate of interest, eq. (2.1), 
to an equilibrium correlator, eq. (5.14) (see ref. [7] for a general argument concerning such 
relations). In order to achieve this goal, the logic is to use an "effective" classical picture to 
describe the long time physics of chemical equilibration. The parameters of this description 
are subsequently matched to reproduce quantum-mechanical correlators. As we will see, the 
consistency of the description will be tested at the matching stage. 

As discussed in sec. 2, large deviations from an equilibrium value tend to decrease, with a 
rate that we want to determine (cf. eq. (2.1)); however, small deviations can also be generated 
by the occasional inverse reactions. 2 This is formally the same physics as in Brownian motion, 
described by a Langevin equation, 

Sh(t) = -r chcm 5n(t)+d(t), (5-11) 

««tU(f)» = n«ten*(t-f). «£(*)» = o, (5.12) 

where 5n is the non-equilibrium excess; £ is a stochastic noise, whose autocorrelation function 
is parametrized by ^^em' an d ((•••)) denotes an average over the noise. The noise is uncorre- 
cted because the time scale considered is much larger than any others in the system. 3 
Now, eq. (5.11) can be solved explicitly, given an initial value Sn(to): 

5n(t) = 5n{t ) e - r <*<»(*-*°) + f* dt > e r chem (t'-t)^/) _ ( 513 ) 

Jto 

2 In a heavy ion collision there may not be enough time for inverse reactions to take place in practice; 
but that does not change the theoretical role that they play in relating the non-equilibrium problem to a 
corresponding equilibrium one. In other words, within the linear response regime the value of the coefficient 
Tchcm is independent of initial conditions and of for how long we observe the dynamics. 

3 At very short time scales, the noise is no longer white but has a structure. By definition, the structure can 
be resolved by inspecting the spectral function corresponding to the "force-force" correlator. As demonstrated 
in sec. 6, the spectral function has support down to small frequencies, with an overall magnitude Q c hem ~ 
e -2M/T Noise becomes coloured at a frequency scale lo vv above which the shape of the spectral function 
changes from its small-frequency asymptotics. This is related to the physics of colour-electric fields, so wc 
may expect uj vv > a 2 g T. This is much larger than the frequency scales that we are concerned with, and plays 
no role in the following. 
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Making use of this solution and taking an average over the noise, we can determine the 2-point 
correlation function of unequal time fluctuations of Sn : 

A cl (t,t') = lim ((Sn(t)Sn(t'))) 

t[)—>— oo 

ft ft' 
= lim / dhe^^-V / dt 2 e r ^ {t2 - t ">{{^{t 1 )^{t 2 ))) 

i-t ft' 

= Q chcm lim / dhe^^-V / dt 2 e r ^ {t2 - t ' ) 5(t 1 -t 2 ) 
to-s-oo J tQ J tQ 

_ ^clicm g-r c hem|t— *'| ^ ]^ 

^-f^chcm 

The limit to — ► — °o here guarantees that any initial transients have died out; therefore, A c i 

is an equilibrium correlation function. Subsequently, making use of dtdf\t — t'\ = —2S(t — t'), 
we obtain 

d t d t ,A cl (t, t') = -^E^l e -r C he m |t-t'l + Qchem s{t _ t ') . (5 . 15) 
Fourier transforming eqs. (5.14) and (5.15) leads to 

/OO Q 
dte^-^A^t') = 2 , (5.16) 

-oo £J + 1 chcm 

/oo 2 o 

dte^(*-*')a t ^A cl (t,t') = ; h 2 cm . (5.17) 

-oo U> ~r L chem 



It is also useful to note that, setting the time arguments equal, we can define a "susceptibility" 

as 

<(5n) 2 ) cl = lim ((6n(t)5n(t))) = , ( 5 .18) 

t -S"OO 21 chem 

where we made use of eq. (5.14). 

Combining eqs. (5.16)— (5.18), various strategies can be envisaged for determining the quan- 
tity that we axe interested, in, namely the non-equilibrium rate r , c j iem . A particularly fruitful 
approach is to take eqs. (5.17), (5.18) as starting points, obtaining 

^chem = r lim ^A cl ( W ), (5.19) 

r c hem < W < <^UV 

rchem = 2((«5n)2)ci ' (5 ' 20) 

Here uj vv is a frequency scale at which some microscopic physics which is not described by the 
effective classical picture sets in, typically o; uv ~ a 2 T, and it has been assumed (cf. sec. 2) 
that r chcm is parametrically small compared with co vv . In our case this is so because r chem 
is exponentially suppressed as ~ e~ M / T . With this input, all real-time information is in the 
numerator of the equilibrium correlator u> 2 A c i(a;). 
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After these preparatory steps, we can promote the determination of r chcm to the quantum 
level. It just remains to note that since in the classical limit observables commute, a suitable 
quantum version of the equilibrium correlator is 

A qm (t,t') = (^{Sh(t),dh(t')}) . (5.21) 

So, eqs. (5.19), (5.20) can be rephrased as 

O chem = Inn u 2 rdte^-^/Ush^Snit')}), (5.22) 

or 

IK - lim rdte^Wi/^) ^)\\ (523) 

chcm " r^S^vl^ 6 \2l dt ' dt' //' (5 - 23) 



together with 



r c hem " 2 {(^7 2 > ' (5 ' M) 



The denominator of eq. (5.24) is nothing but the variance, ((5h) 2 ) = (n 2 ) — (h) 2 . The 
consistency of the matching is tested at least to some extent by whether the variance is 
UV-finite (for most composite operators this is not the case). 

The formulae introduced can be applied on a non-perturbative level by re-expressing them 
through the imaginary-time formalism. This means that we first define a Euclidean correlator, 
fi(r); Fourier-transform it, fl(oj n ) = J^dr e* WnT $7(r), where oj n = 2nnT, n £ % (this requires 
the presence of an UV regulator, or the subtraction of short-distance divergences); and obtain 
the spectral function from its imaginary part, p n {oj) = ImSl(aj ra —> —i[uj + i0 + ]). The sym- 
metric combination needed in eq. (5.23) is given by £\h em = li m r chcm < u < w uv 

The argumentation above can directly be transported to the case at hand, with h replaced 
by H from eq. (5.4). Denoting by Ej the Euclidean electric field, which contains an additional 
i from a Wick rotation, the imaginary-time correlator referred to above reads (we divide by 
volume in order to define intensive quantities) 



n(r) ee i(ftff(T)ftff(o)) v 



= -9 2 [ ( [HE^\ (r, x) [H°E k il>] (0, 0)) , (5.25) 

J x \ I qc 

where gEj, = i[D T , D^], and (...) qc refers to connected quark contractions (the reason for this 
choice is discussed in fig. 3). Hats have been left out in the second row because this correlator 
can be evaluated with regular path integral techniques. Similarly, the correlator related to 
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energy fluctuations becomes 
1 



A(r) ee 



V 



H(t)H(0) 



- IM 



u j + M )ip 



(r,x) 



(0,0) 



(5.26) 



where (...) c refers to the connected part, i.e. (H(t)H(0)) c ee (H(t)H(O)) - (H(0)) 2 . We can 
interpret A(r) as the susceptibility needed in eq. (5.24) to the extent that it is r-independent 
and therefore finite at r — > (cf. eq. (4.1)); this turns out to be the case in the limit irT <C M, 
where it corresponds to a quasi-conserved quantity: A(r) « y((5H) 2 ) (cf. eq. (5.31)). 



5.3. Heavy quark limit 

The correlators in eqs. (5.25), (5.26) can be understood physically, and also written in some- 
what simpler forms, if two-component spinors corresponding to non-relativistic degrees of 
freedom are employed. We choose a representation for the Dirac matrices with 



7°- 



1 




v ~ 1 J 

where are the Pauli matrices. The Dirac spinors are written as 



k = 1,2,3 , 



(5.27) 



IP EE 



X 



4>^{0 ] , -x ] ) 



(5.28) 



Clearly corresponds to and \ to P-^, w hh the projection operators defined as P± = 
\ (l ± 7°). With this notation the operator entering eq. (5.25) can be expressed as 



(5.29) 



Note that this operator is different from that relevant for heavy quark kinetic equilibration: 
electric fields appear in both cases but here they come together with O^Xi X^> whereas in 
ref. [19] the combinations 0^0, x^X appeared. Eq. (5.26) can also be expressed in the new 
notation, with the Hamiltonian becoming 



H 



M{tfe - x ] x) - \ {fa ■ fix + xV • 13 e) 



(5.30) 



For a proper physical interpretation, it is useful to change the ordering of x* a -, Xp- It then 
becomes clear that x* represents an antiparticle to 0; a most direct way to see this is from 
the number density operator: ^7°^ = '4>(P + — P-)ip = 0^0 + X^X = $Q ~ X*^X*- What this 
implies is that operators of the types Q^Xi X ®i appearing in eq. (5.29), create or annihilate 
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quark-antiquark pairs; and that the leading term of the Hamilton operator in eq. (5.30) counts 
particles plus antiparticles, assigning each energies given by their rest mass. 

After these remarks we can simplify the correlator A(r) of eq. (5.26). In the heavy-quark 
limit the leading term comes from M{9^9 — X x) m ec l- (5.30). But since in the same limit 
the cross term gives no contribution, the (disconnect part of) the 2-point correlator is the 
same as that for ?/ry ^ = 9' 9 + x^X- So, 

A(r) « M 2 X{ = M 2 ^((^ 7o ^)(r,x)(V; 7 o^)(0,0)) T , (5.31) 

where X f 1S from eq. (4.1). As required, eq. (5.31) is independent of r. Unfortunately, for 
il(r) of eq. (5.25), it is not clear to us whether any similar simplification is possible; the 
reasons for this are discussed at the beginning of sec. 6. 

To summarize, from the Euclidean correlator, f)(r) in eq. (5.25), we can in principle con- 
struct the Matsubara representation, tl(uj n ) = J dr e* WnT $7(r), if an ultraviolet regulator or 
subtraction is present. After analytic continuation, p n (uj) = lmfl(uj n — > —i[uj + i0 + ]), the 
decay rate of eq. (5.24) follows from 

r = ^Q1^S= lim flhM} (5 32) 

T ^= 2 X{ M- -J^ + W f W' (5 - 32) 

We remark that since eq. (5.25) involves composite operators for non-conserved quanti- 
ties, the issue of renormalization is non-trivial. Unfortunately a satisfactory discussion goes 
beyond the scope of the present work. 



6. Perturbative evaluation 

So far we have made no approximation based on the weak-coupling expansion. At high T, 
however, the renormalized gauge coupling can be assumed small; we would like to make use 
of this limit in order to compare the general formulae with those in sec. 3. 

It is now important to be more precise about the nature of the heavy-quark limit. Even 
though we made use of the "non-relativistic" spinors 9 and X m sec - 5.3 in order to ob- 
tain a physical interpretation for the operators appearing, the function fi(r) cannot actually 
be evaluated with non-relativistic kinematics. A trivial reason is that with non-relativistic 
dispersion relations, a heavy quark and antiquark can annihilate into a single gluon; this 
non-sensical reaction would spoil the physics. In addition, in the t and w-channel processes of 
fig. 1 the heavy quarks are deeply virtual, cf. eq. (3.7). That said, some parts of the analysis 
can still be simplified, but a priori the quark propagators need to be fully relativistic. 

The relevant graphs are shown in fig. 3. It is easy to see that the leading-order graph, (a), 
does not contribute: after analytic continuation and taking the cut we are faced with the 
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decay of a heavy quark and a heavy antiquark into a gluon, which is forbidden by relativistic 
kinematics. At NLO, in contrast, there are non-vanishing contributions; let us show this 
explicitly by evaluating the fermionic graph in fig. 4. 

To get started, we note that in its original form the amplitude squared of eq. (3.4) reads 



Ar fEI-^2l 2 = 9 A N i Tv[T a T b ]Ti[T a T b } 

Tr [7^1 Y?2 ]Tr [7^ + M) 7 „(ff 2 - M)} 



(6.1) 



where T a are the Hermitean generators of SU(-/V C ), normalized as Tr [T a T b ] = whereas 
the imaginary time diagram of fig. 4 can be written as 

n^\cu n ) = -g 4 N f Tr[T a T b ]Tr[T a T b ] 

< 5_K + Pi + P 2 - Ki - K 2 ) e^ a (Pi + P2)e v -APi + P2) 
{PiP2KiK2} PlPKKl + M*) { Kl + M*) 

(Pi + P 2 ) 4 • ^ ' 

Here four-momenta and Dirac-matrices are Euclidean; oj n within the <5 is a short-hand for 
(w n ,0);<5 is normalized so that f± P d(P) = 1; sum-integrals are standard, with $ { ... } denoting 
fermionic Matsubara frequencies; and 

£»;a( P ) = P <V* ~ Pfi $0a (6.3) 

originates from the electric fields. A close kinship between eqs. (6.1), (6.2) is immediately 
observed, but to see that they really lead to the same physics requires a careful analysis. 

We note, first of all, that the index [i appearing in eq. (6.3) can only be spatial. Therefore, 
in the heavy-quark part 

Tr + M) lv (i$ 2 - M)] = 4[V(^i • K 2 - M 2 ) - K^K 2l/ - K lv K 2p \ , (6.4) 

we can drop the terms —K\^K 2v — K\ v K 2li and the spatial part of K\ • K 2 , because the heavy 
quarks will be non-relativistic, cf. eq. (3.5). The part containing final-state momenta, 

{Pi + P2)Su-A p i + ^)Tr [ 7a (if! ) lp (iP 2 )] 
= 4e J;a (P 1 + P 2 )EiAPi + P2) [6 a pPi ■ P2 - Pi a P2/3 - PipP2*\ , (6.5) 

can in turn be re-expressed as (Pj = (p n i, pi)) 

ei;a£i;p6 a p = 3(Pi + P 2 f - 2( Pl + p 2 ) 2 , (6.6) 

£i;a £i;P PlaP2/3 = (Pi + P2) 2 PnlPn2 ~ Pf(Pnl + Pn2)Pn2 ~ P 2 (Pnl + Pn2)Pnl ■ (6.7) 
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(g) 00 (i) 0) (k) (y) 

Figure 3: The graphs contributing to the correlator Q(r) defined in eq. (5.25), up to 0(g 4 ) (time 
runs vertically). The double lines denote heavy quarks; the small dots the composite operators; and 
the grey blob the 1-loop gauge field self-energy. Graphs (a)-(k) look similar to those relevant for 
computing the correlator yielding the heavy quark kinetic equilibration rate [32], but the kinematic 
regime is different. The additional graphs (x) and (y) amount to a renormalization of the gluonic part 
of the energy-momentum tensor by virtual heavy quarks, and have been excluded from the definition 
in eq. (5.25) by restricting to connected quark contractions. 




Figure 4: The part of diagram (k) of fig. 3 sensitive to light quarks, after a Fourier transformation 
to Euclidean frequency u) n and a rotation by 90 degrees. The diagonal line indicates a cut. 



The latter two terms of eq. (6.7) do not contribute due to the antisymmetry in one of the 
summation variables (for instance, in the middle term, after first carrying out T^2 Pnl the 
expression is antisymmetric in p n 2), so we get 



K) « -Sg A C F N c N { ^ d(co n + P 1 + P 2 -K 1 -K 2 )^ 



k nl k n2 - M 2 



1 



f_3 3 2(pi + P2) 2 _ (Pi + P2) 2 2p„ip„ 2 



P?\2P$ (Pi + P2) 2 (Pi + P2Y ~ PWi + P2) 2 ~ Pi (Pi + P2) 2 J ' l ' J 
To carry out the Matsubara sums, we write 

5(u n + Pnl + p n2 - k nl - k n2 ) = / dr e ^+Pm+Pn2-k nl -k n2 )T (6 9) 

Jo 



Then, 



(k nl k n2 - M 2 )e- i ^ +k ^ T 1 e (P-r)(E kl+ E k2 ) +e r(E kl +E k2 ) 

{£jL} ( fc nl + ^)fe + ^ 2 2 ) * 2 (e^ 1+ l)(e^ 2+ l) ' (6 ' ° } 
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where we again approximated E ka « M in the spin part (but not in the exponential functions), 
whereby the "crossed terms" cancelled in the sum. As far as the second row of eq. (6.8) is 
concerned, we note that in the 2nd and 3rd terms a shift p n 2 — > p n 2 — Pni factorizes the 
p ra i-dependence from the r-dependence. These terms lead to a vanishing contribution to the 
transport coefficient defined in eq. (5.32); the reason is that since neither e Pl nor e P2 appears 
in the time dependence, we are left with the phase space constraints S(E kl + E qkl — q) or 
5(Ek 1 + E qkl + q), where E qkl = \/(q — ki) 2 + M 2 and q = pi + p2- These constraints 
cannot get realized and so the factorized terms can be omitted. 4 

Non-trivial contributions arise from the remaining three terms of eq. (6.8). Defining 

rP (P-T){E kl +E k2 ) , r(E k +E k ) p i(p n i+Pn2)r 



^2 (Wn) 



e (P-r)(E kl +E k2 ) _|_ e r(E kl +E k2 ) e i(p nl +p n2 )r ^ + 



rP (p-T)(n kl +n k2 ) , T(H kl +H, k2 ) 

= / dr e itJnT - — T 2 

Jo (e^i +l)(e^ fe 2 +1) r ^ 



+ 1)(^ + 1) {p ^ n2} p i 2p 2 2 ( Pi + p ^ 



(6.12) 

fP . P (P-T)(E k +E k ), T(E k +E k ) J(p n i+Pn2)r 

analytically continuing /9j(u;) = ImXj(u; n — >■ — i[uj + iO + ]); taking the limit uj — > 0; and keeping 
only the terms that give a non-vanishing contribution, some work leads to 

T Pl (u) /p(6 Pi )/ f (6 P2 )[1 - / F (£; fcl )][l - U{E k2 )} 
Jim — — = _ 2ir6(e Pl + e P2 - E kl - E k2 ) , 

(6.14) 

Km _ ljm TpjM x (p,W _ (6 . 15) 

w^o+ cj w^o+ w (pi + p2) z - (e Pl + e P2 ) 

lim = lim x Z!£l^ . (6 . 16 ) 

w-»o+ cj w^o+ w (pi + P2) - (e Pl + e P2 ) z 

In the non-relativistic limit, M 3> 7rT, the subsequent spatial integrals can also be carried 
out. Indeed detailed balance, 

/ F (e pi )/ F (ep 2 )[l - / F (£ fcl )][l - U{E k2 )]5{e Pl + 6 p2 - £7 fcl - E k2 ) 
= f F (E kl )f F (E k2 )[l - / F (e Pl )][l - / P (e P2 )]5(e pi +e P2 - E kl - E k2 ) , (6.17) 

guarantees that the momenta ki,Jt2 are non-relativistic, like in eq. (3.5). Momentum conser- 
vation requires that pi + p2 is also non-relativistic, and that /p(e p J are exponentially small. 



4 In the case with the "double pole" , i.e. the 3rd term of eq. (6.8), one can replace (Pi + P 2 ) 2 — >■ (Pi + P2) 2 + 
m 2 ,; consider first a single pole; and take subsequently a derivative with respect to m\. The relevant phase 
space constraint becomes S(E kl + E qkl — e q ), with e q = ^/q 2 + m\. This does not get realized if mo < 2M, so 
the function vanishes exactly in this regime, and thereby the derivative vanishes as well. 
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So, from eqs. (6.8)-(6.17), 

Iim U^M „ ^ Cf nm I M&Mii&l 

u^o+ oj ^pip 2 kik 2 4e pl e P2 



x (2.) V 4 ) (v, + v 2 - - /c 2 ) { | - 2e ; e ; 2 2 } 

g 4 CWV c iV f r MEki)MEk2) (2vr) 4 ^ 3 )( Pl + p 2 )«5(2 Pl - 2M) 

!/ JpiP 2 kik 2 

/ f F (E kl ) [ U{Ek 2 ) ■ (6.18) 

Jki Jk 2 



2vr 

Dividing by Xf from eq. (4.3), eq. (5.32) finally yields 



ichcm 87TM2 I 27T J 6 ' lb ' iyj 

This agrees with eq. (3.8). 

As far as the gluonic contributions are concerned, the situation is complicated by the many 
diagrams appearing in fig. 3; indeed we have checked that all diagram classes, with two, three 
and four heavy quark propagators, need to be summed together in order to obtain gauge- 
independent results. Nevertheless, without getting lost in excruciating detail, we can draw 
on eqs. (6.1), (6.2) to present a short but "suggestive" argument that things work out as 
before. For the s-channel process, the vacuum amplitude squared reads 

= g ^Tv[T a T b }f acd f bcd If (?i)P«?(P 2 ) 

TrirWi + M) 7 "(# 2 -M)} 
x {V l +V 2 Y 

X [Vap{V 2 ~ Vl)p ~ T] plx {Vl + 2V 2 )a + VnaPPl + V 2 )p\ 

x [vap{V 2 - Vi) v - Vpu{Vi + 2V 2 ) d + TfrtiZPi + V 2 )p] . (6.20) 

Here Pr denotes the projector from a sum over the on-shell gluon polarizations, and Feynman 
gauge was used for the inner gluon line. On the other hand, the gluonic equivalent of the 
process in fig. 4 can be written in Feynman gauge as 

59, {99) {u n ) = A g ^Ti[T a T b ]f acd f bcd 

gK + P 1 + P 2 -K 1 - K 2 ) e WQ (P 1 + P 2 )e^(Pi + P 2 ) 
PiP2{KiK2} P?P?(K! + M*)(K1 + Mi) 

Trfr^ff! +M) lv (i$2 — M)\ 
(Pi + P 2 ) 4 
x [S a p(iP 2 - iPi) a &pa 

(iPi + 2iP 2 ) a + S aa (2iP 1 + iP 2 ) p ] 
x [5 ap {iP 2 - iPi)p - StfiiPi. + 2iP 2 ) a + 8p a {2iP 1 + iP 2 ) p ] . (6.21) 
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Establishing a precise equivalence between all indices requires adding other gluonic contri- 
butions on both sides, but a comparison with eqs. (6.1), (6.2), for which we carried out a 
detailed analysis, allows us to anticipate that things work out here as well, including the 
important factor \ in front of the gluonic channels in eq. (3.3), clearly visible in eq. (6.21). 



7. Discussion 

The question of whether or not heavy quarks chemically equilibrate in heavy ion collisions is 
sometimes addressed by comparing the observed total yield with that predicted by a thermal 
distribution at the final (pionic) freeze-out temperature. In this paper, we have have asked 
whether chemical equilibrium could be reached earlier on, at a higher temperature. Since 
there are many heavy quarks in the initial state, one simply needs to get rid of some of them, 
to arrive at a thermal ensemble. The rate for this is suppressed by e~ M//T , which is the 
density of antiquarks seen by any given heavy quark. If this suppression can be overcome 
then, for a while, heavy quarks could be part of the thermal medium, before re-decoupling 
again above the final pionic freeze-out, explaining why more heavy quarks and antiquarks are 
observed than is due for chemical equilibrium. 

Taking the expression from eq. (3.10); factorizing from it the susceptibility of eq. (4.3); 
normalizing the susceptibility to its value in the massless limit, to be denoted by Xo = N c T 3 /3; 
and setting N c = 3, the result for the chemical equilibration rate reads 

Setting furthermore iVf = 3, a s ~ 0.3, M ~ 1.5 GeV, and estimating Xf/Xo from refs. [27, 28], 
we obtain T" 1 ^ ~ 10 fm/c at T ~ 600 MeV, and r^ m > 60 fm/c at T ~ 400 MeV. If true, 
these time scales indicate that chemical equilibrium is unlikely to be reached in current heavy 
ion collision experiments, where the highest temperatures are around T ~ 400 MeV and the 
time scale is around 10 fm/c. 

The estimate presented in eq. (7.1) is a rough one. In principle, a non-perturbative value 
could be obtained from eq. (5.32) through numerical lattice Monte Carlo simulations and 
a subsequent analytic continuation. For the latter step, short-distance singularities need 
to be subtracted, as has recently been elaborated upon in connection with other transport 
coefficients [33, 34]. This task is undoubtedly a hard one: as an analysis of graph (a) of fig. 3 
shows, for oj » M the spectral function behaves as 

Pn{uj) = I20(4vr)3 ^ + ° ( - UJ M >J ' (7 ' 2) 

implying that the Euclidean correlator diverges as f2(r) ~ 1/t 7 for r <C M~ 1 . To subtract 
this dominant and any subdominant divergences perturbatively, and still retain a statistically 
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significant signal containing the thermal physics, would require a very precise analysis. (Al- 
ternatively one could start with the correlator A(r) of eq. (5.26), although this is dominated 
by a constant mode, which poses problems for some methods of analytic continuation.) 

Nevertheless, our non-perturbative formulation may have other uses; for instance, it may 
be amenable to an order-of-magnitude estimate in the confined phase through chiral effective 
theories, similarly to what has previously been achieved in the case of the heavy flavour kinetic 
equilibration rate [35]-[38]. Possibly it could also be combined with non-relativistic QCD 
(NRQCD) where the hard (p ~ M) momentum fields have been integrated out perturbatively. 
Indeed it is possible to include the effects of QQ annihilation in NRQCD, through a 4-fermion 
interaction in the effective Lagrangian, where the effective coupling has an imaginary part [39]. 
In this case one cannot consider fi(r) of eq. (5.25) because the chromo-electric field is hard 
and should have been integrated out; but one could compute A(r) of eq. (5.26) instead. 

We end by remarking that whereas our non-perturbative formulation is only valid near 
equilibrium, the Boltzmann description can also be applied beyond it. Since r chem is propor- 
tional to the density of the antiquarks, cf. eqs. (3.1)-(3.3), we may expect a correspondingly 
faster rate in the real world where the heavy antiquarks appear in overabundance. 
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